---
title: "Week 8 Homework"
author: "David Barron"
date: "Trinity Term 2018"
output: pdf_document
---

```{r, echo=FALSE, message=FALSE, warning=FALSE}
knitr::opts_knit$set(width = 100)
knitr::opts_chunk$set(comment =  '', echo=FALSE)
knitr::opts_chunk$set(fig.align = 'center', fig.width = 6, fig.asp = 0.618, out.width = '80%')
options('show.signif.stars'=FALSE)
options('digits' = 3)
#options(scipen=5)
library(ggplot2)
library(foreign)
library(lavaan)
library(semPlot)
#library(psych)
```

```{r funcs}
my.mi <- function(model){
  mi <- lavaan::modificationindices(model)
  ix <- order(-mi[, 4])
  mi <- mi[ix, ]
  ix <- mi[, 4] > 1
  if (sum(ix) <= 10) return(mi[ix, ])
  else return(mi[1:10, ])
}
```
  
```{r homework}
ch10.cov <- matrix(c(1.024, 0, 0, 0, 0, 0, 0, 0, 0, 
                     0.792, 1.077, 0, 0, 0, 0, 0, 0, 0,
                     1.027, 0.919, 1.844, 0, 0, 0, 0, 0, 0,
                     0.756, 0.697, 1.244, 1.286, 0, 0, 0, 0, 0, 
                     0.567, 0.537, 0.876, 0.632, 0.852,  0, 0, 0, 0,
                     0.445, 0.424, 0.677, 0.526, 0.518, 0.670, 0, 0, 0,
                     0.434, 0.389, 0.635, 0.498, 0.475, 0.545, 0.716, 0, 0,
                     0.580, 0.564, 0.893, 0.716, 0.546, 0.422, 0.373, 0.851, 0, 
                     0.491, 0.499, 0.888, 0.646, 0.508, 0.389, 0.339, 0.629, 0.871),
                   ncol = 9, byrow = TRUE)

colnames(ch10.cov) <- c('EdAsp', 'OcAsp', 'VerbAch', 'QuantAch', 'FamInc', 'FaEd', 'MoEd', 'VerbAb', 'QuantAb')

print(ch10.cov)

```         
The model is that a student's achievement (`Achieve`) depends on home and family characteristics (`Home`), the student's ability (`Ability`), and the student's aspiration (`Aspire`).  These four latent variable have the following measures:

_Achieve_: Verbal achievement (VerbAch) and Quantitative achievement (QuantAch).
_Home_: Family income (FamInc), Father's education (FaEd) and Mother's education (MoEd).
_Ability_: Verbal ability (VerbAb) and Quantitative ability (QuantAb).
_Aspire_: Educational aspiration (EdAsp) and Occupational aspiration (OcAsp)

The covariance matrix above is based on data from 200 students.

Aspiration is thought to be influenced by Ability and Home situation.  Achievement is thought to be influenced by all three latent variables.  Specify an appropriate structural equation model, obtain results and determine whether the fit of the model is adequate.  If it is not, modify the model and re-fit.

  
  
# First SEM

```{r mod}
ch10.mod1 <- ' Home =~ FamInc + FaEd + MoEd
Ability =~ VerbAb + QuantAb
Aspire =~ EdAsp + OcAsp
Achieve =~ VerbAch + QuantAch

Aspire ~ Home + Ability
Achieve ~ Home + Ability + Aspire

Home ~~ Ability
'
ch10.op1 <- sem(ch10.mod1, sample.cov = ch10.cov, sample.nobs = 200)
summary(ch10.op1)
fitMeasures(ch10.op1, c('chisq', 'df', 'pvalue', 'rmsea', 'cfi', 'srmr'))
semPaths(ch10.op1, what = "est")

my.mi(ch10.op1)
residuals(ch10.op1, type = 'standardized')$cov
```
There are a number of possible modifications that these results suggest.  Examples include covariance between father's and mother's education (which seems very plausible).  Let's try this out.

```{r model}

ch10.mod2 <- ' Home =~ FamInc + FaEd + MoEd
Ability =~ VerbAb + QuantAb
Aspire =~ EdAsp + OcAsp
Achieve =~ VerbAch + QuantAch

Aspire ~ Home + Ability
Achieve ~ Home + Ability + Aspire

Home ~~ Ability

FaEd ~~ MoEd
'

ch10.op2 <- sem(ch10.mod2, sample.cov = ch10.cov, sample.nobs = 200)
summary(ch10.op2)
fitMeasures(ch10.op2, c('chisq', 'pvalue', 'rmsea', 'cfi', 'srmr'))

semPaths(ch10.op2, what = "paths")
```
That seems to do the trick.  You can see that all the hypothesised structural relationships are statistically significant.